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ANALYSIS  OF  LAYERED  COMPOSITE  PLATES  BY 
THREE-DIMENSIONAL  ELASTICITY  THEORY 


Abstrjiqt 

This  report  contains  two  parts:  the  first  part  is  concerned  with 
the  natural  vibrations  of  laminated  anisotropic  plates,  and  the  second 
part  is  concerned  with  geometrically  nonlinear,  finite-element  analysis 
of  the  bending  of  cross-ply  laminated  anisotropic  composite  plates. 
Individual  laminae  are  assumed  to  be  homogeneous,  orthotropic  and 
linearly  elastic.  A  fully  three-dimensional  isoparametric  finite  element 
with  eight  nodes  (i . e . ,  linear  element)  and  24  degrees  of  freedom  (three 
displacement  components  per  node)  is  used  to  model  the  laminated  plate. 
The  finite  element  results  of  the  linear  analysis  are  found  to  agree  very 
well  with  the  exact  solutions  of  cross-ply  laminated  rectangular  plates 
under  sinusoidal  loading.  The  finite  element  results  of  the  three- 
dimensional,  geometricall^  nonlinear  analysis  are  compared  with  those 
obtained  by  using  a  shear  deformable,  geometrically  nonlinear,  plate 
theory.  It  is  found  that  the  deflections  predicted  by  the  shear  defor- 
mable  plate  theory  are  in  fair  agreement  with  those  predicted  by  3-D 
elasticity  theory;  however,  stresses  are  found  to  be  not  in  good  agree¬ 
ment.  The  results  of  natural  vibration  indicate  that  for  relatively 
thick  plates,  the  shear  deformable  plate  theory  predicts  frequencies 
higher  than  those  predicted  by  the  3-D  theory. 
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NATURAL  VIBRATIONS  OF  LAMINATED  ANISOTROPIC  PLATES 
USING  3D-ELASTICITY  THEORY 
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and 

T.  Kuppusamy 
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ABSTRACT 

The  paper  contains  a  description  of  the  three-dimensional 
elasticity  equations  and  the  associated  finite  element  model  for  natural 
vibrations  of  laminated  rectangular  plates.  The  numerical  results  for 
natural  frequencies  are  compared  with  those  obtained  by  a  shear 
deformable  plate  theory.  A  number  of  cross-ply  and  angle-ply 
rectangular  plates  are  analyzed  for  natural  frequencies.  For  relatively 
thick  plates,  the  plate  element  predicts  frequencies  higher  than  those 
predicted  by  the  3-D  element. 


INTRODUCTION 


The  present  study  is  motivated  by  the  lack  of  three-dimensional 

finite-element  results  for  natural  frequencies  of  laminated  composite 

plates.  For  relatively  thick  plates  (l.e.,  for  side-to-thlckness  ratios 

less  than  10),  the  classical  plate  theory  and  certain  thick  plate 

theories  predict  higher  frequencies,  and  therefore  It  Is  of  Interest  to 

have  a  3-D  finite  element  to  accurately  predict  natural  frequencies. 

The  following  literature  provides  the  background  for  the  present  study. 

The  shear  deformable  plate  theory  (see  Reissner  [1]  and  Mindlin 

[2])  is  a  two-dimensional  theory  derived  from  the  assumption  that  the 

strains  e  ,  e  ,  and  e  are  independent  of  z,  which  leads  to  the 
z  xz  yz 

displacement  field  of  the  form 

u  «  uQ(x,y)  +  z<j»x (x,y) 
v  *  vQ(x,y)  +  z*y{x,y) 

w  =  wQ(x,y).  (1) 

Here  u0,  v0,  and  w0  denote  the  displacements  along  x,  y,  and  z 

directions  respectively  of  a  point  located  on  the  midplane  of  the  plate, 

which  coincides  with  the  xy-plane  of  the  coordinate 

system;  4 >  and  4<  denote  the  rotations  (taken  clockwise  positive)  of  a 
x  y 

line  element,  initially  perpendicular  to  the  midplane,  about  y  and  x 

axes,  respectively.  The  classical  plate  theory  is  a  special  case 

(derived  from  the  assumption  that  e  =e=e=0)of  the  shear 

z  xz  y  z 

deformable  plate  theory  in  which  the  transverse  shear  strains  are 


assumed  to  be  zero: 
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The  classical  plate  theory  Is  based  on  the  Ki rchof f-Love  assumption  that 
planes  Initially  normal  to  the  midplane  remain  plane  and  normal  to  the 
midsurface  after  bending,  l.e.,  equivalently,  conditions  In  Eq.  (2) 
hold. 

The  first  lamination  thin  plate  theory  Is  due  to  Relssner  and 
Stavsky  [3].  An  extension  of  the  thick  plate  theory  to  arbitrarily 
laminated  anisotropic  plates  was  presented  by  Yang,  Norris,  and  Stavsky 
[4].  Whitney  and  Pagano  [5],  and  Reddy  and  Chao  [6]  presented  exact 
solutions  of  the  theory  when  applied  to  laminated  rectangular  plates 
under  certain  lamination  schemes,  boundary  conditions,  and  sinusoidal 
distribution  of  transverse  loads.  Reddy  and  his  colleagues  [7-12] 
presented  f Ini be -element  analyses  of  the  bending,  vibration  and 
transient  response  of  laminated  anisotropic  composite  plates. 

The  three'-dlmenslonal  lamination  theory  Is  based  on  the  assumption 
that  the  individual  lamina  behave  according  to  the  3-D  elasticity 
theory.  At  the  Interfaces  of  Individual  lamina,  the  displacements  are 
assumed  to  be  continuous.  Dana  and  Barker  [13]  used  the  cubic 
isoparametric  brick  element  and  Putcha  and  Reddy  [14]  used  the  linear 
brick  element  to  investigate  the  bending  of  layered  plates.  Recently, 
Kuppusamy  and  Reddy  [15]  used  the  linear  Isoparametric  brick  element  to 
study  the  geometrically  nonlinear  behavior  of  laminated  plates.  The 
present  study  is  concerned  with  the  application  of  the  linear  brick 
element  to  the  natural  vibration  of  laminated  composite  plates. 

GOVERNING  EQUATIONS 

Consider  a  plate  arbitrarily  laminated  of  a  finite  nuitber  of 
orthotropic  layers  of  uniform  thickness.  The  equations  of  motion  of  a 
three-dlmensltonal  elastic  continuum  under  the  assumption  of 


Infinitesimal  strains.  In  the  Lagrangian  description,  are  given  by 


.00 


3  boj  . 

z  — I* 


j=i  aXj. 


*2n(k) 

(k)  a  ul  , 

P  — 2 —  »  1 

at 


1,2,3, 


(3) 


where  <?jk)  *  denote  the  components  of  the  Cauchy  stress 
tensor,  Is  the  density,  and  u|k^  are  the  displacement  components  of 
the  k-th  layer.  All  of  the  variables  are  referred  to  the  plate 
coordinates  x^.  It  Is  assumed  that  there  are  no  body  forces  in  the 
problem.  At  the  interface  of  two  layers  It  is  assumed  that  the 
displacements  are  continuous. 


uOO  „  U(k+1) 

U1  ui 

The  strain-displacement  relations  are  given  by 


(4) 


5U. 


I  /  I  A 

eij  “  2  'dXj 


au . 


ax. 


i»j  =  1,2,3. 


(5) 


The  stress-strain  relations,  for  Individual  lamina,  in  the  plate 
coordinates  are  given  by 


,(k)  = 


6 

Z 

j  =  l 


c(k>ejk> 
ij  j 


1  *  1,2, 


,6, 


(6) 


where 

a,  =  on,  a2  *  a22,  o3  =  O33,  o4  =  0,.  =  a^,  ag  - 

®1  *  ell»  e2  *  e22’  e3  *  e33’  e4  *  2e23*  e5  *  2e13*  e6  *  2e12’(7) 

and  cW  are  the  material  stiffnesses  of  k-th  layer  referred  to  the 
plate  coordinates. 
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Equations  (2),  (3),  and  (6)  are  to  be  solved  using  appropriate 
boundary  conditions.  At  a  point  on  the  boundary  of  the  plate,  one 
should  specify: 

A 

either  u^  =  u^  (1  =  1,2,3), 

or  nj °i j  s  (1  =  1,2,3),  (8) 

A 

where  nj  denotes  the  j-th  component  the  unit  outward  normal  n  to  the 

A  A 

boundary,  and  u^  and  t^  denote  specified  values  of  the  displacement  and 
traction  components. 


FINITE-ELEMENT  MODEL 

The  variational  formulation  of  Eqs.  (2),  (3),  and  (6),  over  a 
typical  element  Qe,  is  given  by 


a2u. 


0  * 


I#  *  p6ui  dxdy'12  ■  VijV5 


(9) 


where  6  denotes  the  variational  symbol.  Let  u^  be  interpolated  over  the 
element  Qe  by 
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u.  =  E  uU  (1  =  1,2,3)  (10) 

1  a=l  '  a 

where  $  (x,y,z)  denote  the  trilinear  Interpolation  functions  of  the 
a 

eight-node  isoparametric  brick  element,  and  u*(t)  denote  the  value  of  u^ 
at  node  a  and  time  t.  Substitution  of  Eq.  (10)  into  Eq.  (9),  we  obtain 


[M11]  [M12]  [M13] 

..  \ 
l«!> 

[K11]  [K12]  [K13] 

[M22]  [M23] 

i 

{u2} 

>+ 

[K22]  [K23] 

< 

{ u2^  > 

33 

„  33, 

symm.  [M  J 

.M 

symm.  IK  ] 

{ U31 
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where  [K^j]  and  [M^j]  are  coefficient  matrices  associated  with 
vector  {u.}  in  the  i-th  equation.  For  the  sake  of  brevity  the  explicit 
form  of  the  coefficients  is  not  given  here.  The  mass  and  bending 
coefficients  are  to  be  evaluated  using  full  integration  and  the  shear 
coefficients  (in  [K])  are  to  be  evaluated  using  the  reduced  integration. 

To  investigate  the  free  vibration,  we  set  { F^. }  *  { 0} ,  and 

u?  =  %ePWt  (i  =  1,2,3;  a  *  1,2 . 8)  (12) 

where  p  =  /T  and  id  denotes  the  frequency  of  natural  vibration. 
Substitution  of  Eq.  (12)  into  Eq.  (11)  gives  the  following  generalized 
eigenvalue  problem: 

([K]  -  u>2[M]){A}  =  {0},  (13) 

where  (a}  =  U<J10}.  { ,  {u30}}T.  Equations  (13)  are  assembled  in  the 
usual  manner,  and  (homogeneous)  boundary  conditions  are  imposed  before 
solving  the  eigenvalue  problem  by  any  standard  eigenvalue  program. 

NUMERICAL  RESULTS 

In  all  of  the  examples  considered  here,  the  biaxial  symmetry  is 
exploited  to  model  only  one  quadrant  of  the  plate.  The  quadrant  is 
modeled  using  4  x  4  x  n  mesh  of  the  eight-node  brick  elements.  Here  n 
denotes  the  number  of  layers  in  the  plate.  The  following  two  types  of 
high-modulus  composite  materials  are  used  in  the  present  study: 

Material  I:  E^E^S,  E^E^l,  <WGi3/y0-5’  G23/E2*0,2 

V12*V13’V23‘°*ZS*  l,“-0 
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Material  II:  Ej/Eg-40,  E3/E2=l,  G12/VS13/E2=0*6,  G23/E2=0.2 

v12=v13=v23=0.25,  p=1.0.  (14) 

.  o 

A  value  of  5/6  is  used  for  the  shear  correction  coefficient,  iq.  The 

boundary  conditions  used  in  the  present  study  are  shown  in  Fig.  I. 

The  first  example  Is  concerned  with  three-layer  cross-ply  square 
plates:  0°/90°/0°,  h^  =  h3  =  h/4,  h2  =  h/2.  This  is  also  equivalent  to 
four-layer  cross-ply  (0°/90o/90o/0°)  laminate  made  of  equal  thickness 
layers.  Table  1  contains  a  comparison  of  the  fundamental  frequencies 
obtained  by  the  plate  theory  [6]  and  the  present  3-D  elasticity 
theory.  It  is  clear  from  the  results  that  the  solutions  predicted  by 
the  shear  deformable  plate  theory  are  higher  than  those  predicted  by  the 
3-Delasticity  theory;  of  course,  the  classical  thin  plate  theory 
predicts  even  higher  frequencies.  The  difference  can  be  explained  as 
follows:  in  plate  theories,  the  transverse  shear  strains  and  normal 
strains  are  either  completely  neglected  (like  in  the  classical  plate 
theory)  or  included  in  an  approximate  sense  (e.g.,  the  transverse  shear 
strains  are  included  in  the  shear  deformable  plate  theory). 

Consequently,  the  modeled  plate  is  stlffer  than  the  actual  one.  Due  to 
the  lower  value  of  the  shear  modulus  relative  to  the  in-plane  Young's 
moduli,  the  transverse  shear  deformation  effects  are  more  pronounced  in 
composite  plates.  In  the  3-D  theory,  no  assumption  is  made  to  neglect 
the  shear  or  normal  strains,  and  therefore  the  frequencies  predicted  are 
more  accurate  (and  realistic).  For  large  ratios  of  side  to  thickness, 
the  transverse  shear  strains  and  normal  strains  are  negligible  and 
therefore  both  theories  predict  almost  the  same  frequencies.  It  is  also 
clear  from  the  results  presented  in  Table  1  that  the  degree  of  material 


r 

7 


(a)  Boundary  conditions  for  cross-ply  plates(BCl) 


(b) 


Boundary  conditions  for  angle-ply  plates  (BC2) 


Figure  1  Boundary  conditions  used  In  the  present  study  for 
cross-ply  and  angle-ply  plates 
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anisotropy  as  well  as  plate  aspect  ratio  adversely  affect  the  accuracy 
of  the  frequencies  predicted  by  the  plate  theory. 

The  second  example  Is  concerned  with  the  vibration  of  a  square 
sandwich  plate.  The  material  of  the  face  sheets  Is  the  same  as  Material 
I,  and  that  of  the  core  Is  transverselylsotropic  with  the  following 
properties: 


E:/E2  =  1.0,  E3/E2  =  12.5,  p  =  1.0 

G13  =  G23  =  1,5E2*  G12  =  0,4E2 
v12  =  0.25,  v13  =  v23  =  0.02.  (15) 

The  thickness  of  each  face  sheet  is  0.1 h  and  that  of  the  core  Is  0.8h. 


Table  2  contains  the  nondimensionalized  frequencies  obtained  using  the 
3-D  elasticity  theory.  These  results  should  serve  as  bench  marks  for 


future  comparisons. 

Table  3  contains  nondimensionalized  natural  frequencies  of  two- 
layer  cross-ply  (0°/90°)  and  angle-ply  (45°/-45°)  square  plates.  From 
the  results  presented  In  Table  3  and  4,  it  follows  that  two-layer  cross- 
ply  plates  have  lower  frequencies  than  both  two-layer  angle-ply  plates 
and  four-layer  cross-ply  (0o/90o/90o/0°)  plates  (of  equal  thickness 
layers).  On  the  other  hand,  the  two-layer  angle-ply  plates  have  higher 
frequencies  than  the  cross-ply  plates  analyzed  here.  These  observations 
indicate  that  the  two-layer  cross-ply  laminates  are  structurally  more 
flexible  whereas  the  two-layer  angle-ply  laminates  are  more  stiffer  than 
the  four-layer  (0o/90°/90o/0°)  cross-ply  plates. 

The  results  of  the  natural  vibration  of  four  different  laminations 


of  four-layer  angle-ply  square  plates  are  discussed  next.  Table  4 
contains  the  nondimensionalized  frequencies  of  the  following  four  cases 


of  laminates: 
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Table  1.  Nondimen  si  onal  1  zed  fundamental  frequency,  \  =  /p/Eg,  of 
three-layer  cross-ply  (0°/90°/0°),  simply  supported  plates 


(hi  3  h3  =  h/4,  h2  =  h/2). 


Material 

I 

Material 

II 

Finite  Element 

Method 

percentage 

Finite  Element  Method 

percentage 

b/a 

a/h 

Plate  Theory  3-0  Theory 

error 

Plate  Theory 

3-D  Theory 

error 

5 

8.831 

8.317 

-6.18 

9.896 

9.119 

-8.52 

1 

10 

12.380 

11.805 

-4.87 

14.451 

13.370 

-8.09 

(1 2.233)t 

(14 .295)  t 

100 

15.474 

15.473 

-0.0065 

19.177 

18.959 

-1.15 

3 

10 

- 

- 

- 

12.555 

11.338 

-10.734 

tobtained  using  2x2  mesh  of  nine-node  elements. 


Table  2.  Nondimensional ized  frequencies  of  a  square  sandwich  plate  with  simply 
supported  boundary  conditions,  (b/a  =  1 ;  h^  *  h3  =  O.lh,  h2  =  0.8h; 
30-Theory 


Mode 

a/h 


10  11.132 


2  3  4  5  6  7 


32.881  37.153  71.787  78.006  79.616  83.380 


Table  3.  Nondlmensionallzed  fundamental  frequencies  of  two-layer  (0°/90°  or 
-45°/45°)  square  plates  with  simply  supported  boundary 
conditions  (a/b  -  1,  a/h  -  10,  hj  -  h/2;  30-Theory) 


Lamination 

Scheme 

1 

Material  I 

2  3 

4 

1 

Material  II 

2  3 

4 

0°/90<’ 

9.176 

44.094  44.044 

60.886 

10.563 

48.806  48.806 

67.032 

45°/-45° 

16.668 

33.816  36.723 

46.929 

17.826 

37.690  39.260 

51.305 

Table  4.  Nondimenslonal ized  fundamental  frequencies  of  angle-ply  square  plates 
(material  I,  a/h  »  10) 


_ -45V_45°/-45° 

3D  Theory 

Mode 

3D 

Plate 

45°/0°/45° 

-30o/30o/-30° 

-30°/0°/-30° 

1 

16.483 
(1 6.546)t 

15.282 

18.246 

16.602 

16.141 

2 

34.044 

(36.129) 

34.100 

36.565 

34.513 

33.615 

3 

35.051 

(46.273) 

45.062 

46.484 

45.975 

48.984 

4 

45.753 

(64.354) 

46.591 

49.281 

47.845 

52.255 

ft he  numbers  in  parenthesis  correspond  to  the  case  in  which  the  inplane 
Inertias  are  omitted. 


1.  (-45 °/45 °/-45° ) ,  hj  =  h3  =  h/4,  h?  =  h/2,  a/h  *  10 

2.  (-45°/0°/-45°).  hj  =  h3  =  h/4,  *  h/2,  a/h  =  10 

3.  (-30°/30°/-30°),  h1  =  h3  =  h/4,  h?  =  h/2,  a/h  =  10 

4.  (-30°/0°/-30°) ,  hj  =  h3  =  h/4,  h2  =  h/2,  a/h  =  10 

The  material  properties  used  are  those  of  material  I.  Table  4  also 
contains  natural  frequencies  for  Case  1  when  the  inplane  Inertias  are 
omitted.  It  is  clear  that  the  effect  of  the  inplane  inertia  is  to 
reduce  the  frequencies,  and  this  reduction  has^signi ficant  effect  on 
higher  modes. 
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A  THREE-OIMENSIONAL  NONLINEAR  ANALYSIS  OF  CROSS-PLY 
RECTANGULAR  COMPOSITE  PLATES 
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Virginia  Polytechnic  Institute  and  State  University 
Blacksburg,  VA  24061 

Abstract 

The  results  of  a  three-dimensional,  geometrically  nonlinear, 
finite-element  analysis  of  the  bending  of  cross-ply  laminated 
anisotropic  composite  plates  are  presented.  Individual  laminae  are 
assumed  to  be  homogeneous,  orthotropic  and  linearly  elastic.  A  fully 
three-dimensional  isoparametric  finite  element  with  eight  nodes  (i.e., 
linear  element)  and  24  degrees  of  freedom  (three  displacement  components 
per  node)  is  used  to  model  the  laminated  plate.  The  finite  element 
results  of  the  linear  analysis  are  found  to  agree  very  well  with  the 
exact  solutions  of  cross-ply  laminated  rectangular  plates  under 
sinusoidal  loading.  The  finite  element  results  of  the  three- 
dimensional,  geometrical ly  nonlinear  analysis  are  compared  with  those 
obtained  by  using  a  shear  deformable,  geometrically  nonlinear,  plate 
theory.  It  Is  found  that  the  deflections  predicted  by  the  shear 
deformable  plate  theory  are  in  fair  agreement  with  those  predicted  by  3- 
n  elasticity  theory. 
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1.  INTRODUCTION 

Composite  materials  exhibit  higher  stiffness-to-weight  ratios  and 
increased  corrosion  resistance  compared  to  isotropic  materials.  The 
anisotropic  material  properties  of  layered  composites  can  be  varied  by 
varying  the  fiber  orientation  and  stacking  sequence.  While  this  feature 
gives  the  designer  an  added  degree  of  flexibility,  the  stiffness 
mismatch  of  the  orthotropic  layers  bonded  together  with  different  fiber 
orientations  leads  to  interlaminar  stresses  in  the  vicinity  of  free 
edges.  Tor  certain  stacking  sequences,  loading,  and  boundary  conditions 
interlaminar  stresses  can  be  so  large  that  they  dictate  the  design  of 
the  structure. 

Analyses  of  layered  composite  plates  can  be  divided  into  two 
groups:  (i)  analyses  based  on  a  laminate  plate  theory,  and  (ii) 
analyses  based  on  a  three-dimensional  laminated  (elasticity)  theory. 

The  laminated  plate  theory  Is  the  extension  of  the  classical  plate 
theory  (CPT)  or  the  Rei ssner-Mindl in  shear-deformable  plate  theory  (SDT) 
to  layered  composite  plates.  The  first  lamination  theory  including 
hending-stretching  coupling  is  apparently  due  to  Reissner  and  Stavsky 
[11.  Yang,  Norris,  and  Stavsky  [?]  presented  a  generalization  of  the 
Reissner-Mindl in  thick-plate  theory  for  homogeneous,  isotropic  plates  to 
arbitrarily  laminated  anisotropic  plates.  Whitney  and  Pagano  [31  (also 
see  Reddy  and  Chao  [4])  presented  closed-form  solutions  to  the  theory 
when  applied  to  certain  cross-ply  and  angle-ply  rectangular  plates. 

Reddy  [5]  presented  a  finite-element  analysis  of  the  lamination 
theory.  A  higher-order  lamination  theory  that  accounts  for  a  cubic 
variation  (as  opposed  to  linear  in  [2-5])  of  the  inplane  displacements 
and  quadratic  variation  of  the  transverse  displacement  through  the 
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thickness  was  presented  by  Los  Christensen,  and  Wu  [6],  and  hybrid- 
stress  finite-element  analysis  of  the  theory  was  presented  by  Spilker 
[7]. 

In  laminated  plate  theory  it  is  assuned  that  the  laminate  is  in  a 
state  of  plane  stress  (as  assumption  carried  from  the  classical  plate 
theory)  and  integrals  through  the  thickness  of  a  laminate  are  equal  to 
the  sum  of  integrals  through  the  thickness  of  individual  laminae.  These 
assumptions  lead  to  inaccurate  prediction  of  interlaminar  stresses  at 
the  free  edges,  although  the  solution  is  reasonably  accurate  away  from 
free  edges.  The  laminate  plate  theory  is  not  accurate  in  a  boundary 
layer  region  which  extends  inward  from  the  edge  to  a  distance 
approximately  equal  to  the  laminate  thickness. 

The  fully  three-dimensional  laminate  theory  Is  an  extension  of  the 
elasticity  theory  of  a  three-dimensional  solid  composed  of  layers  of 
different  material  properties.  Pipes  and  Pagano  [8]  and  Pipes  [9]  used 
a  finite-difference  technique  to  solve  the  quasi-three-dimensional, 
linear,  elasticity  equations  for  laminates  (also  see  Hsu  and  Herakovich 
[10]).  Lin  [11],  Oana  [12],  and  Oana  and  Barker  [13]  used  a  cubic, 
three-dimensional,  isoparametric  element  with  72  degrees  of  freedom  to 
analyze  laminated  plates  (also  see  Putcha  and  Reddy  [14]).  The 
numerical  results  in  these  studies  agree  very  well  with  those  of  Pagano 
[15,18]. 

The  present  study  Is  motivated  by  the  lack  of  finite-element 
results  for  three-dimensional,  geometrically  nonlinear  analysis  of 
layered  anisotropic  composite  plates.  A  finite-element  formulation  of 
the  geometrically  nonlinear  theory  of  a  laminated,  three-dimensional, 
elastic  contlnuun  is  presented.  Numerical  results  of  the  linear  as  well 
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as  nonlinear  analysis  are  presented  for  several  cross-ply  plate 
problems.  The  formulation  is  validated  by  comparing  the  present  results 
of  the  linear  analysis  with  those  of  Pagano  [15,16].  The  results  of  the 
nonlinear  analysis  should  serve  as  bench  mark  results  for  future 
investigations. 


2.  GOVERNING  EQUATIONS 

Consider  a  laminate  (C)  composed  of  N  orthotropic  layers  with  axes 
of  elastic  symmetry  parallel  to  the  plate  axes.  The  laminate  is 

A 

subjected  to  normal  traction  t^  =  q(x^,x,,)  at  its  upper  surface  (i.e., 
X3  =  h/2).  The  constitutive  equations  for  any  layer  are 
given  by 
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where  'oij-  and  Ty  are  the  components  of  the  stress  and  strain  tensors, 
respectively,  defined  in  the  material-coordinates,  and  TT^.  are  the 
material  stiffness  coefficients.  The  coefficients,  ,  are  given 

in  terms  of  the  engineering  constants  by  (see  Jones  [17]) 

^11  =  V1  "  V23V32^/A 
ri2  *  E2^v12  +  v32v13*A 
^13  =  E3*v13  +  v1Pv23^/A 
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C22  *  M1  "  v13v31^/A 
^23  =  E3^v23  +  v21v13^/A 
C33  =  E3^  "  v12v21^/A 
^44  =  fi23  *  ^55  =  G13  »  ^66  =  G12 
A  =  (]  "  v12v21  '  v23v32  ‘  v31v13  '  ?v21v31v32^’ 

where 


Ej.E2.E3  *  Young's  nodul i  in  1,  2  and  3  (material)  directions, 
respectively 

Vj  .  =  Poisson's  ratio  for  transverse  strain  in  the  J- 

direction  when  stressed  in  the  i-direction 

G23,G13*^12  =  Shear  moduli  in  the  2-3,  1-3,  and  1-2  planes, 

respectively.  (3) 


In  view  of  the  reciprocal  relations 


(4) 


there  are  only  nine  independent  elastic  constants  for  an  orthotropic 
elastic  medium. 

The  constitutive  equations  (1)  when  transformed  to  the  plate 
coordinate  system  take  the  form 
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where  C^j  are  the  transformed  plate  stiffnesses. 

The  nonlinear  kinematic  description  of  an  elastic  body  yields  the 
following  equations  of  equilibria  (in  the  absence  of  body  forces  and 
moments) 


5u 


J  1 


(6) 


wherein  denotes  the  Kronecker  delta  symbol,  and  the  summation 


e  =  1  fill  +  .^1  +  122  ZE) 
u  ?.  lax,  ax,  ax,  ax,;* 


(7) 


convention  on  repeated  subscripts  is  used.  The  strain-displacement 
equations  of  the  large-deflection  theory  of  elasticity  are  given  by 

Du,  DUj 

■i  ~‘t  *  5xi  ’‘"3 
To  complete  the  description  of  the  field  equations,  Eqs.  (4)-(6), 
and  (7)  should  he  adjoined  by  boundary  conditions.  At  any  point  of  the 
boundary  of  the  body  one  should  specify  one  of  the  following  two  types 
of  boundary  conditions: 

(i)  essential  (or  geometric)  boundary  conditions 

A 

u_  =  u 
m  m 

(ii)  natural  (or  dynamic)  boundary  conditions 

t  =  n,<j. ,  (6.  +  um  ,)  =  tm.  (8) 

m  j  1  j '  mi  m,i '  m  '  ' 

Here  nj  denotes  the  j-th  component  of  the  unit  normal  to  the  boundary, 

A  A 

tm  the  m-th  component  of  the  boundary  traction,  and  um  and  tm  denote 
specified  values. 


3.  FINITE-ELEMENT  FORMULATION 
Here  we  present  a  displacement  finite -element  model  of  the 
equations  (3)-(8).  To  this  end  we  construct  the  variational  formulation 
of  the  equations  over  an  arbitrary  element  of  the  finite-element 
mesh.  We  have  (see  [18,  p.  382]) 
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Here  6  denotes  the  variational  symbol,  and  is  given  in  terms  of 
via  Eqs.  (4)  and  (7). 

The  displacements  are  interpolated  by  expressions  of  the  form. 


8  „ 

um  =  2  (m  =  1»2>3)  (10) 

a=l 

where  <|>  (x,y,z)  (a  =  1,2,.. .,8)  are  the  trilinear  interpolation 
oc 

functions  of  the  eight-node  isoparametric  element  in  three  dimensions, 
and  u“  denotes  the  value  of  um  at  node  a.  Substituting  Eq.  (10)  into 
(9) ,  we  obtain 
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F?  -  /  t  *  ds  ,  (j,m,n  =  1,2,3)  (12) 

“  r(e)  m  8 

Every  isoparametric  finite  element  of  the  finite-element  mesh 
can  be  generated  from  the  master  element  via  the  transformation  (see 
Figure  2) 

xi  3  X1  ^1*52,53^ 
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Z  xi^a(^l »^2 ’^3^  ^  = 

o=l 


(13) 


where  x“  are  the  global  coordinates  of  the  element  nodes 

and  (i  =  1,2,3)  are  the  local  coordinates.  Therefore,  the  integrals 

in  Eg.  (12)  can  be  transformed  to  the  master  element  and  evaluated 


numerically  using  the  Gauss  quadrature.  The  transformation  of 
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For  example,  consider 


9<Ji  54-fi  11  1  .  . 

J  (e)  'SxT'SxT  dxldx2dx3  =  (^l»^^3^1dli2dC3 
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where  Pj  and  Wj  are  the  Gauss  points  and  weights,  respectively ,  and  the 

Integrand  F1^  (i,j  =  1,2,3;  a,p  =  1,2,... ,8)  is  given  by 
op 
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(16) 


Kk£  being  the  elements  of  the  inverse  of  the  Jacobian  matrix,  [J], 

This  procedure  can  be  implemented  on  a  digital  computer,  and  the  element 
coefficient  matrices  in  Fq.  (10)  can  be  evaluated  numerically.  It  is 
well  known  that  when  the  ratio  of  side  to  thickness  of  the  plate  is  very 
large  (i.e.,  when  the  plate  is  thin,  say  a/h  >  20)  one  should  use 
reduced  integration  to  evaluate  the  shear  terms  (i.e.,  terms  involving 
C44,  C45,  and  C55).  For  a  trilinear  element,  Ixlxl  Gauss  quadrature 
must  be  used  to  evaluate  the  coefficients  of  O44,  O55,  and  2x2x2  Gauss 
rule  to  evaluate  all  other  terms. 

Since  the  coefficient  matrix  [K]  depends  on  the  unknown  solution 
vector  {u},  one  should  employ  an  iterative  solution  procedure  to  solve 
the  finite-element  equations.  Here  we  use  the  Picard  type  iterative 
technique,  which  begins  with  an  assumed  displacement  field  (usually,  set 
to  zero  to  obtain  the  linear  solution)  to  compute  [K]  at  the  beginning 
of  the  first  iteration.  In  subsequent  iterations,  solution  obtained 
from  the  previous  iterations  is  used  to  compute  [K],  The  iteration  is 
continued  until  the  solutions  obtained  in  two  consecutive  iterations 
differ  by  a  preassiqned  error  margin  (say,  10-4).  It  is  more  economical 
to  use  load  incremental  methods  in  conjunction  with  the  iterative 
technique  described  above.  In  other  words,  for  each  increment  of  the 
load  the  increment  to  the  nonlinear  solution  is  obtained  (see  [IP]). 
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4.  DISCUSSION  OF  THE  NUMERICAL  RESULTS 
Results  of  the  Linear  Analysis 

In  order  to  validate  the  present  formulation  and  element,  first  the 

linear  analysis  of  a  symmetric  three-layer  square  laminate  with 

the  7-^-di rection  (of  material  principal  axes)  coinciding  with  the  x^ 

direction  (of  the  plate  axes)  in  the  outer  layers  and 

the  X2~direction  parallel  to  the  xi-axis  in  the  center  layer  is 

performed.  The  thickness  of  the  outer  layers  is  assumed  to  be  one-half 

of  the  thickness  of  the  center  layer  (h1  =  hj  =  h/4,  h2  =  h/2).  The 

loading  is  assumed  to  he  sinusoidal  with  respect  to  the  x-|-x9  plane, 

*  TtX.  -Itx  p 

t3  5  <l(xi»x2>  =  q0  cos  ~a~  cos  ~b~  (17) 


and  the  boundary  conditions  are  of  the  simply  supported  type  which  allow 
normal  displacement  on  the  boundary,  but  prevent  tangential 
displacement.  For  a  quarter  plate  these  imply 
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for  the  engineering  constants 


Et  =  1.724  x  10R  kN/m2  (25  x  10fi  psi) 
E?  *  E3  =  fi.89  x  ID6  kN/n2  (1(1*  psi) 
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G12  =  G13  =  3*4S  x  1o6  kN/m2  (°*S  x  1()6  Ps1)  (19) 

G23  =  1  3.78  x  106  kN/m2  (0.2  x  106  psi) 


v12  =  V31  =  v32  =  n**5 


The  following  nondimenslonal ization  Is  used  to  present  the 
displacements  and  stresses: 
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foj  *  °2,0i2^  =  2  ^al  ,CJ2*a12^  ’  ^23,a13^  =  q  a  ^23,013^ 
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The  nondimensionalized  center  deflection  H3(0, 0,0)  obtained  using  the 
laminated  plate  theory  [4,5]  and  the  elasticity  theory  are  presented  in 
Table  1.  One  can  conclude  from  the  results  that  the  finite  element 
results  obtained  by  using  the  reduced  integration  (R)  are  in  good 
agreement  with  the  analytical  solution  of  Pagano  [16],  and  the  shear- 
deformable  plate  theory  solution  [4,5]  differs  from  the  3-0  elasticity 
solution  by  about  1 0%  for  the  problem  at  hand.  The  classical  plate 
theory  (OPT)  solution  differs  from  the  shear  deformable  plate  theory 
(SOPT)  by  35%  for  side  to  thickness  ratio  of  a/h  =  10. 

The  nondimensionalized  stresses  (a^,  "cig  =  0^3.  05  =  o^2)  for  the 
problem  are  compared  in  Table  2.  We  observe  that  for  a/h  >  50  the 
results  obtained  using  the  reduced  integration  are  in  good  agreement 
with  the  analytical  solutions,  and  for  a/h  <  20  the  full  integration 
gives  better  results. 

Nondimensionalized  center  deflection  and  stresses  of  three-layer 
(h^  =  h/3)  rectangular  plates  (b/a  =  3)  of  cross-ply  (0°/90o/0°) 
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Table  1  Nondimensionalized  center  deflection^  in  a  square  cross-ply 
(0°/90°/0°)  plate  (h,  =  h3  3  h/4,  h2  =  h/2)  under  sinusoidal 
loading  (Eg  =  E^ j  Ei  =  2 5 & ^ *  ^12  *  ^13  =  ^*23  s  ^*^2* 

v12  '  v23  '  v13  ■  °-25)  "  ‘  " 


Shear 

Elasticity  Solution  Deformable  Plate  Solution[4] 

a/h  Pagano  2x2x3  4x4x3  2x2x3  4x4x3  CFS  2x2  4x4  2x2  4x4  CPT 
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0.431 

100 

0.435 

0. 

.027 

0.089 

0, 

.412 

0.429 

0.434 

0.032 

0 

.103 

0, 

.412 

0.428 

0.431 

+  u3(0,0,»);  mesh  shown  is  for  a  quarter  plate;  boundary  conditions  are 
tne  same  as  those  shown  in  Eq.  (18). 
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Table  2.  Comparison  of  nondimenslonal  stresses  *  for  three-layer 

(0°/90°/0°)  cross-ply  square  plate  under  sinusoidal  loading 


(a/b  =  1 

1.  hj  = 

hg  =  h/4,  1*2  = 

h/2) . 

a 

¥ 

30-El asti city  Theory 
Stress  Pagano  2x2x3  4x4x3  2x2x3 
type  [16]  F  F  R 

4x4x3 

R 

Plate  Theory  [4] 
2x2  4x4  2x2 

F  F  R 

4x4 

R 

*1 

1.388 

0.603 

0.676 

0.133 

0.165 

0.259 

0.318 

0.259 

0.318 

2 

°5 

0.153 

0.224 

0.224 

0.179 

0.197 

0.242 

0.286 

0.255 

0.289 

"e 

0.086 

0,047 

0.053 

0.025 

0.031 

0.027 

0.033 

0.028 

0.033 

0.720 

0.389 

0.456 

0.185 

0.232 

0.302 

0.377 

0.306 

0.379 

4 

°5 

0.219 

0.182 

0.209 

0.138 

0.205 

0.276 

0.330 

0.293 

0.335 

°6 

0.047 

0.029 

0.034 

0.021 

0.025 

0.023 

0.029 

0.024 

0.029 

*1 

0.599 

0.326 

0.425 

0.268 

0.331 

0.334 

0.451 

0.380 

0.467 

10 

°5 

0.301 

0.235 

0.281 

0.147 

0.286 

0.316 

0.389 

0.353 

0.400 

0.028 

0.028 

0.022 

0.015 

0.018 

0.016 

0.022 

0.019 

0.023 

‘l 

0.543 

0.240 

0.391 

0.298 

0.364 

0.265 

0.437 

0.404 

0.494 

20 

*5 

0.328 

0.421 

0.342 

0.277 

0.312 

0.303 

0.395 

0.373 

0.420 

®6 

0.023 

0.010 

0.017 

0.013 

0.016 

0.011 

0.018 

0.017 

0.021 

*1 

0.539 

0.083 

0.234 

0.308 

0.376 

0.098 

0.279 

0.413 

0.504 

50 

*5 

0.337 

0.755 

0.990 

0.286 

0.321 

0.252 

0.348 

0.380 

0.427 

°6 

0.021 

0.003 

0.009 

0.012 

0.015 

0.004 

0.011 

0.016 

0.020 

*1 

0.271 

0.025 

0.095 

0.310 

0.378 

0.030 

0.120 

0.414 

0.505 

100 

"5 

0.339 

0.877 

1.569 

0.287 

0.274 

0.230 

0.298 

0.381 

0.428 

0.339 

0.001 

0.004 

0.012 

0.015 

0.001 

0.005 

0.016 

0.020 

*The  stresses  In  the  finite-element  analysis  are  computed  at  the  Gauss 
points. 
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Table  3.  Nondlmenslonallzed  center  deflection  In  a  rectangular  (b  =  3a) 
cross-ply  (0°/90°/0°)  plate  (h,  =  h,  =  h,  =  h/3)  under  sinu¬ 
soidal  loading  (BC1 ,  E2  =  E3,  E,  =  25E?,jG19  -  G13  =  0.5E,, 

G23  -  °-2E2>  »„  *  -  0.35)  2  12  13  2 

-  - Shfflf - 

Elasticity  Solution  Deformable  Plate  Solution  [5]  CPT 

a/h  Pagana  2x2x3  4x4x3  2x2x3  4x4x3  CFS  2x2  4x4  2x2  4x4  Solution 


2  8.17  7.953 

4  2.82  2.702 

10  0.919  0.801 
20  0.610  0.411 
50  0.508  0.137 
100  0.503  0.042 


7.831  8.410  7.934 
2.716  2.905  2.764 
0.872  0.915  0.902 
0.540  0.589  0.600 
0.305  0.494  0.513 
0.134  0.481  0.501 


10.11  8.153  7.900 
2.97  2.418  2.376 
0.93  0.748  0.788 
0.60  0.425  0.530 
0.51  0.150  0.327 
0.50  0.050  0.153 


8.600  8.000  0.480 
2.655  2.405  0.480 
0.810  0.804  0.480 
0.560  0.574  0.480 
0.490  0.509  0.480 
0.480  0.500  0.480 
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Table  4.  Comparison  of  nondlmenslonalized  stresses*  for  three-layer 
cross-ply  (0°/90°/0°)  square  plate  under  sinusoidal  loading 
(a/b  =  1,  ht  -  h/3) 


a 

h 

Stress  Pagano 
Type  [1 5] 

30-El astl city  Theory 
2x2x3  4x4x3  2x2x3 

F  F  R 

4x4x3 

R 

2x2 

F 

Plate  Theory  [5] 

4x4  2x2  4x4 

F  R  R 

*1 

0.938 

0.550 

0.619 

0.097 

0.122 

0.274 

0.337 

0.273 

0.336 

2 

0.309 

0.257 

0.267 

0.212 

0.236 

0.124 

0.295 

0.263 

0.290 

°6 

0.070 

0.044 

0.050 

0.025 

0.032 

0.037 

0.046 

0.038 

0.046 

®1 

0.755 

0.400 

0.470 

0.162 

0.205 

0.324 

0.406 

0.330 

0.408 

4 

°5 

0.282 

0.221 

0.247 

0.209 

0.244 

0.284 

0.340 

0.303 

0.345 

*6 

0.051 

0.031 

0.036 

0.021 

0.025 

0.028 

0.034 

0.030 

0.035 

0.590 

0.325 

0.424 

0.244 

0.301 

0.343 

0.464 

0.392 

0.481 

10 

*5 

0.357 

0.271 

0.333 

0.301 

0.343 

0.311 

0.382 

0.348 

0.393 

°6 

0.029 

0.017 

0.022 

0.014 

0.017 

0.017 

0.023 

0.020 

0.024 

0.552 

0.233 

0.381 

0.267 

0.327 

0.267 

0.440 

0.408 

0.499 

20 

*5 

0.385 

0.426 

0.361 

0.327 

0.369 

0.295 

0.381 

0.360 

0.405 

0.023 

0.010 

0.016 

0.012 

0.014 

0.011 

0.019 

0.017 

0.021 

*1 

0.541 

0.080 

0.225 

0.274 

0.335 

0.098 

0.280 

0.414 

0.504 

50 

°5 

0.393 

0.757 

0.992 

0.335 

0.377 

0.248 

0.338 

0.363 

0.408 

°6 

0.022 

0.003 

0.009 

0.011 

0.013 

0.004 

0.011 

0.017 

0.020 

®1 

0.539 

0.024 

0.092 

0.276 

0.336 

0.030 

0.120 

0.414 

0.505 

100 

0.395 

0.878 

1.570 

0.337 

0.378 

0.229 

0.290 

0.364 

0.409 

®6 

0.021 

0.001 

0.004 

0.011 

0.013 

0.001 

0.005 

0.016 

0.020 

*The  stresses  In  the  finite-element  analysis  are  computed  at  the  Gauss 
points. 
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construction  are  presented  in  Tables  3  and  4.  The  loading,  boundary 
conditions,  and  material  properties  are  the  same  as  those  given  in  Eqs. 
(17),  (18),  and  (19),  respectively.  The  finite-element  results  are  in 
good  agreement  with  the  corresponding  exact  solutions  [15]  and  plate 
theory  solutions  [5], 

Results  of  the  Non  1  inear  Analysis 

First,  geometrically  nonlinear  analysis  of  isotropic  plates  is 
performed  and  the  results  are  compared  with  the  results  available  in  the 
literature.  A  square  isotropic  plate  (v  =  0.3)  subjected  to  uniformly 
distributed  load  is  analyzed  using  simply  supported  boundary 
conditions.  Due  to  the  biaxial  symmetry,  only  one  quadrant  of  the  plate 
is  modeled  with  4x4x1  mesh  of  linear  elements.  Since  membrane  and 
bending  contributions  dominate  the  stiffness  matrix,  full  integration  is 
used  to  evaluate  the  element  matrices.  The  present  results  for  the 
center  deflection  and  stresses  are  compared  In  Figs.  1  and  2  with  the 
nonlinear  thin  plate  theory  (i.e.,  the  von  Karman  plate  theory) 
solutions  available  in  the  literature  [20,21].  From  the  results  of  the 
isotropic  plate  one  can  conclude  that  the  nonlinearity  exhibited  by  the 
3-0  elasticity  theory  Is  relatively  smaller  than  that  included  in  the 
classical  von  Karman  plate  theory  [21]  but  larger  than  that  in  the  shear 
deformable  plate  theory  [201. 

Next,  a  symmetric  cross-ply  (0°/90°/0°)  square  plate  under 
uniformly  distributed  load  Is  analyzed.  The  boundary  conditions  and 
material  properties  are  given  by  Eqs.  (18)  and  (19),  respectively.  The 
3-0  elasticity  results  for  center  deflection  and  stresses  are  compared 
with  the  shear  deformable  plate  theory  results  [20]  in  Figs.  3-5.  The 
effect  of  thickness  (i.e.,  ratio  of  side  to  thickness)  on  the 
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:i] 


Figure  1  Load-deflection  curves  for  simply  supported  Isotropic 
(v  -  0.3)  square  plate  under  uniform  load. 
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Load  parameter 


'°r2 


Figure  2  Center  normal  stress  versus  the  load  parameter  for 
simply  supported  Isotropic  (v  ■  0.3)  square  plate 
under  uniform  load  (a  Is  computed  at  x  *  (0.0264, 
0.0264,  0.0789  In  3-D*elast1c1ty  theory  and  x  • 
(0.0625,  Q.Q625,  0.1)  In  the  plate  theory) 


•  K  ' 


i 


Figure  3  Load-deflection  curves  for  cross-ply  (0°/90°/0°, 
equal  thickness  layers)  square  plates  under 
uniformly  distributed  load  (E^/E?  *  25,  G,~  * 
S13  *  °-5E2,G23  ■  0.2  E2,  v12'-  0.25. 
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Figure  4  Normal  stress  versus  load  parameter  for 
three-layer  cross-ply  (0°/90o/0o)  square 
plate  under  uniform  loading  (S  *  a/h  *  10). 
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Figure  5  Center  normal  stress  versus  load  parameter 
for  three-layer  cross-ply  (0°/90°/0°)  square 
plate  under  uniform  load  (S  s  a/h  «  4). 


Figure  6  Load-deflection  curves  for  two-layer,  cross- 
ply  (0°/90°)  square  plate  under  uniformly 
distributed  load  (Ei  -  25E?,  G-i?  *  Gn  »  0.5E?, 
G23  -  0.2EZ,  v12  *  0.25). 


3-0  Elasticity 
SDPT 

3-0  Elasticity 


Load  parameter,  P  = 


Figure  7  Normal  and  shear  stresses  versus  load  parameter 
for  two-layer,  cross-ply  (0°/90O,S*10)  square 
plate  under  uniformly  distributed  load  (o  is 
computed  at  x  »  (.0625,  .0625,  0.1)  in  SDPT; 
ov_  is  computed  at  x  a  (.4736,  .4736,  .0894)  in 
30zand  at  x  *  (.4375,  .4375,  0.1)  in  SDPT). 
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deflections  and  stresses  can  also  be  seen  from  Figs.  3-5.  It  should  be 
noted  that  the  plate  theory  solutions  deviate  more  from  the  3-D 
elasticity  solutions  as  the  thickness  increases.  Also,  it  should  be 
pointed  out  that  the  stresses  in  the  two  theories  are  computed  at 
different  locations  of  the  plate,  and  therefore  part  of  the  difference 
in  the  two  solutions  is  attributable  to  the  difference  in  the  location 
of  points  at  which  the  stresses  are  calculated. 

Finally,  a  two-layer  cross-ply  (0°/90°)  square  plate  under 
uniformly  distributed  load  is  analyzed  and  the  plots  of  center 
deflection  and  stresses  versus  the  load  parameter  are  presented  in  Figs. 
6  and  7.  Note  that  the  nonlinearity  exhibited  by  the  3-D  elasticity 
theory  is  less  than  that  of  the  plate  theory  for  load  parameter  values 
less  than  250. 


5.  CONCLUSIONS 

A  finite-element  analysis  of  geometrically  nonlinear,  three- 
dimensional  theory  of  laminated  plates  is  presented.  It  is  found  that 
the  shear  deformable  plate  theory  results  are  fairly  accurate  when 
compared  to  the  three-dimensional  elasticity  theory  for  side-to- 
thickness  ratios  of  10  or  more.  The  difference  between  the  two 
solutions  is  larger  for  a  side-to-thickness  ratio  of  4.  It  is  also 
found  that  the  stresses  predicted  by  the  plate  theory  are  in  larger 
error  than  the  deflections  when  compared  to  those  predicted  by  the 
three-dimensional  elasticity  theory.  The  reduced  integration  technique 
is  found  to  have  an  effect  on  the  accuracy  of  the  solution:  reduced 
integration  is  recommended  for  thin  plates  (a/h  >  10)  and  full 
Integration  for  thick  plates  (a/h  <  10),  especially  when  geometric 
nonlinearities  are  included. 
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The  first  part  of  the  report  is  concerned  with  the  natural  vibrations  of 
laminated  plates,  and  the  second  part  is  concerned  with  geometrically 
nonlinear  analysis  of  cross-ply  laminates  plates.  A  fully  three-dimensional 
linear  isoparametric  element  is  used  to  determine  natural  frequencies,  bending 
deflections  and  stresses. 
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